Vortex lattices in Bose-Einstein condensates: 
from the Thomas-Fermi to the lowest Landau level regime 
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We consider a periodic vortex lattice in a rotating Bose-Einstein condensed gas, where the cen- 
trifugal potential is exactly compensated by the external harmonic trap. By introducing a gauge 
transformation which makes the Hamiltonian periodic, we solve numerically the 2D Gross-Pitaevskii 
equation finding the exact mean field ground state. In particular, we explore the crossover between 
the Thomas-Fermi regime, holding for large values of the coupling constant, and the lowest Landau 
level limit, corresponding to the weakly interacting case. Explicit results are given for the equation 
of state, the vortex core size, as well as the elastic shear modulus, which is crucial for the calculation 
of the Tkachenko frequencies. 
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In the last years, a significant effort has been devoted 
to the study of vortex lattices in harmonically trapped 
rotating Bose-Einstein condensates. Striking results have 
been obtained experimentally, leading to the observation 
of large vortex lattices in fast rotating condensates f 
to the measurement of their dynamical properties 
Theoretically, several predictions have been provided jj, 
0,0 m agreement with experiments and more challenging 
regimes, related to the quantum Hall effect, have also 
been proposed [?|]. 

In this paper we study the rotating analogue of a uni- 
form Bose-Einstein condensate, where a periodic vortex 
lattice is present. For dilute atomic gases, which are 
highly compressible, this can only be obtained when the 
centrifugal potential is exactly cancelled by the harmonic 
confinement, i.e., when the harmonic trapping frequency 
is precisely equal to the angular velocity fi. Under this 
condition, in spite of the non-periodicity of the Gross- 
Pitaevskii (GP) equation governing the system, the den- 
sity n (see Fig.^l and the velocity field v' in the rotating 
frame are periodic in the plane of rotation || , as can 
be proved by a proper gauge transformation. Analyt- 
ical predictions are available in the literature for both 
the weakly interacting regime 0, 0, 0, 0, 0] , where 
the system is in the lowest Landau level (LLL), and the 
strongly interacting one 0] , where the density is basically 
constant and the Thomas-Fermi (TF) approximation ap- 
plies. Solving numerically the GP equation, we recover 
these limiting cases and compute the ground state for 
any value of the interaction strength. The validity of the 
mean-field approach requires that the number of atoms 
per vortex be very large. If this condition is not satisfied 
the system enters a strongly correlated regime related to 
the quantum Hall effect y|- 

In the rotating frame, the external confinement V cx t 
combines with the centrifugal potential, giving rise to the 
effective trapping V c s — V cx t — mfl 2 (x 2 + y 2 )/2, where 
m is the atomic mass. As anticipated above, we impose 
the compensation V c s — and decouple the motion in 



FIG. 1: Density distribution of the condensate in the x-y 
plane computed by solving Eq. @ for (a) g = 0, (b) g = 300. 
In (b) the dashed hexagon is the Wigner-Seitz cell, while the 
solid parallelogram is the computation box (see text). Darker 
regions correspond to lower density. 



the axial direction from the radial one, therefore reduc- 
ing to a 2-dimensional problem with an effective cou- 
pling constant g2D 01- I n addition, our discussion can 
be conveniently reformulated in terms of dimensionlcss 
quantities by using harmonic oscillator units, i.e., Ml, SI, 
and Iq = \JhjrMl for energy, frequency, and length re- 
spectively. The length Iq also fixes the equilibrium inter- 
vortex distance, according to the Feynman expression for 
the vortex density n v — rMl/nh = 1/ttIq |l5| . Having 
to deal with an infinite periodic vortex lattice, we fix 
the average particle density (n) in terms of the number 
iVccii of atoms per lattice cell, (n) = N cc u/A cc u, where 
A;eU = = 7T is the cell area. The corresponding di- 
mensionless stationary Gross-Pitaevskii equation in the 
rotating frame is hence 



(—iV - e z A r) 



(i) 



where /i is the chemical potential, e z is the unit vec- 
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tor along z, and the dimensionless coupling constant 
g = ng2D (n) = g2DN ce nm/h 2 is the only parameter 
governing the equation. The normalization is J \ip\ 2 = 1, 
where J = J. dr means integration over a single cell. 

Due to the presence of the rotational term e z A r, 
Eq. Q is not spatially periodic. However, in the s piri t 
of the magnetic field analogy discussed in Refs. 0, UJ| , 
one can perform a gauge transformation which turns the 
symmetric effective vector potential A = e z A r into a 
periodic function A' = A — VA. We choose A = St, 
where St is the (non-periodic) phase associated with the 
Tkachenko (non-periodic) velocity field vt = VSV used 
by Tkachenko to describe the vortex lattice of an incom- 
pressible fluid [9j. With such a choice, indeed, one finds 
A' = e z A r — Vt = —v' T , where the Tkachenko velocity 
v' T in the rotating frame is now periodic [jj. We there- 
fore write ip = i>e lST and substitute this expression into 
Eq. l[T]l. finally obtaining 



flip 



-iV 



(2) 



which admits periodic solutions. The explicit expres- 
sion of vt is conveniently written using the complex 
variable representation z — x + iy, v — v x + iv y |l6| . 
The lattice geometry is determined by the half-periods 
w\ and u>2 (see Fig. Cfb)) which, due to the constraint 
^4-ccii = I", are related by Im^Ju^) = 7r/4. One then 
has i>t(z*) = «[C( Z ) where ((z) = £(z; u)\, o^) = 

1 /z + y^l- [1 /(z — Zjfe) + 1 /zifc + z/z?,„1 is the quasi-periodic 
Weierstrass Q- function [l7j. Here the primed summation 
excludes the term j — k — 0, Zjfc = 2jwi+2fca;2 are the vor- 
tex positions, and a = [u>* — £(ui)]/uji is chosen to guar- 
antee the periodicity of v' T {z, z*) = vt(z*) — iz. Corre- 
spondingly, the Tkachenko phase is given by St(z, z*) = 
arg a(z) + Im(o!Z 2 /2), where the Weierstrass tr-function 
obeys d z a{z) = ((z)a(z). 

In the following, we will determine the numerical solu- 
tion of the GP equation in the periodic Tkachenko gauge, 
Eq. (J2J, for arbitrary values of g. Since the ground state 
solution of Eq. @ can acquire a non- vanishing phase, the 
velocity field associated with -0 does not coincide in gen- 
eral with Vt, although the deviations are always found 
to be very small. Actually, the velocity field coincides 
with vt not only in the incompressible regime (g — > oo) 
but also in the LLL limit (g — > 0). Indeed the LLL wave 
function can be written in the form 0] 

^lll = ]J(z z jfe ) e -' 2 ' 2 / 2 = ( T(z)e- 2 / 2 e-l^ 2 / 2 (3) 
jk 

and hence, by using the explicit expression for St re- 
ported above, one remarkably gets ^lll = |'0LLL|e 4 ' ST . 

We have solved Eq. by restricting the calculation 
to the parallelogram given by the lattice vectors (see 
Fig. IHb)) and imposing periodic boundary conditions. 
After discretizing the continuous problem on a triangu- 
lar mesh, we used the finite element method |l8| to find 



40 
30 
=L 20 
10 


a, 1.15 

ro 

i 1.1 

ro 

K 1.05 
1 





TF-^ 




LLL 






I ■ 




■ 


J T - LLL 











20 



40 



60 



80 



100 



FIG. 2: Upper panel: numerically computed chemical po- 
tential /i (solid line) as a function of g. The dashed and 
dash-dotted lines are the LLL and TF predictions respec- 
tively. Lower panel: derivative ndfi/dg. Numerical results 
(points with error bars) reproduce the correct behaviour in 
both the LLL and TF regimes. 



solutions by means of imaginary time propagation. The 
energies of different lattices are calculated by varying the 
lattice vectors at constant cell area. We verified that the 
triangular geometry (where a — 0) is the favourite one 
for any value of g > 0, while for g = any lattice with 
^4ccii = I" has the same energy, reflecting the infinite LLL 
degeneracy. For g < we find that the system is un- 
stable, as evident from its negative compressibility. In 
Fig. ^ the equilibrium density profile extracted from the 
numerical solution of Eq. (0) for g — > + and g = 300 is 
presented. Vortices are characterized by holes in the den- 
sity, whose size increases by decreasing the interaction. 
When j > 1, the density is uniform everywhere except 
in a very narrow region of size much smaller than the 
inter- vortex distance and the TF approximation applies. 
On the contrary, if g <§; 1 the vortex core radius becomes 
comparable to Iq and the density profile is highly struc- 
tured (density maxima and saddle points alternate on the 
boundary of the hexagonal Wigner-Seitz cell). 

It is also interesting to study the behaviour of the 
chemical potential /i as a function of g. Making a pertur- 
bative expansion around the LLL triangular lattice, one 
can prove that \i is linear in g for g — > 0, i.e., \i ~ l+0q/ir 
with (5 = /|Vlll|7(J>lll| 2 ) 2 = 1.1596 0. 
By using our numerical solution, the equation of state 
/i = fx(g) can be calculated for any value of g, providing 
the bridge between the g = result and the TF relation 
/i = g/n, valid for g 1. Since these two opposite limits 
both exhibit a linear dependence with a different slope, 
the transition region shows a non linear behaviour. Devi- 
ations from linearity are however very small. The results 
are reported in Fig. [21 where we also plotted the deriva- 
tive dfi/dg to highlight the transition. 
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FIG. 3: Fractional vortex core area A as function of 7J j \ jL . 
Solid line: numerical solution of Eq. Dashed line: LLL 
( f lll -* °°) limit from E q- ©• Dash-dotted line: TF 
(^lll ~* 0) limit A — 4 £ 2 . Points with error bars: experi- 
mental results from Ref. |2^| (see text). 




FIG. 4: Shear modulus C2/(n) as a function of Points 
connected by the solid line: numerical solution of Eq. JS). 
Dashed line: LLL prediction of Ref. fH\ . Thin solid line: TF 
prediction of Ref. |2ql . In the inset the same data are plotted 
as a function of g (bottom label) and T~^ LL (upper label). 



The behaviour of fj,(g) (g oc (n)) can be used to evalu- 
ate the coarse grained density profile (n) of the trapped 
rotating gas in terms of the residual effective potential 
V c ff, through the local density relationship n((n)) + V c g — 
const. In particular, for the experimentally relevant case 
of harmonic trapping, one notices that the quasi-linear 
behaviour of /i(g) implies the density profile to be an 
inverted parabola |2o| . 

Another important quantity to study is the vortex core 
size. In the TF limit this length scale is fixed by the heal- 
ing length £ = l/y/2jl, while in the LLL limit, where £ 
diverges, the core size saturates to a fraction of the inter- 
vortex distance 0, 0] . Several definitions of the vortex 
core radius r v are possible. For a detailed comparison 
with experiments, we follow the procedure of Ref. [22| . 
where r v is defined as the mean square root radius of the 
Gaussian which better fits the function n max — n(r), eval- 
uated on a single lattice cell 23]. Our numerical results 
are reported in Fig. |3 where we plot the fractional core 
area A = irr1/A cc \\ — r% as a function of the inverse of the 
LLL parameter r LLL = (/x- l)/2 22]. For r LLL > 1 the 
vortex core radius follows the TF behaviour (r v ~ 1.98£), 
while for Tlll — > it reaches the limiting value A = 0.34. 
The latter result is not far from the value A = 0.30 deriv- 
able using the Gaussian fit definition with the circular cell 
approximation of Ref. [f| . 

Concerning the comparison with experimental results 
[22|. two data sets are presently available. They differ in 
the procedure employed for the expansion, which is re- 
quired to reach a sufficient resolution in the imaging pro- 
cess. The first set (open squares in Fig.|3J) gives the frac- 
tional core area as measured during an expansion which 
preserves the vertical size, while the second one (filled cir- 



cles) is obtained allowing also for some axial expansion. 
Being similar to a 2D scaling, the first procedure is ex- 
pected to preserve the initial value of the fractional area 
A [3] and indeed the corresponding experimental data 
are in good agreement with theory. The second tech- 
nique, involving an additional axial expansion, brings in- 
stead A closer to the LLL limiting value, since it lowers 
the value of the effective 2D coupling constant. 

We finally studied the elastic shear modulus C%. This 
quantity, related to the energy variation due to a shear 
distortion of the lattice |25J, is a crucial ingredient for 
the determination of the frequencies of Tkachenko modes 
[HIM HIE! El- In order t0 extract this coefficient, we 
calculate the energy for different lattice shapes, starting 
from the triangular lattice and then varying one of the 
lattice vectors at fixed cell area: lu' 2 — > lu2 + ewi. The 
coefficient C2 is then given by C2 = 35£/4e 2 , where S£ 
is the (second order) variation of the energy per lattice 
cell. 

Our numerical results are reported in Fig. 0] as a func- 
tion of Tlll ■ In the TF limit Tlll ,g>lwe recover the 
result Cij{n) = 1/8 [2s|. while in the LLL limit we find 
perfect agreement with the prediction of Refs. [ill l27j. 
Cij(n) ~ 0.1191(7, which is a factor 10 larger than the 
value previously estimated in Ref. |2sj . The figure shows 
that the transition between the two regimes is character- 
ized by a maximum at Tlll — 5. A consequence of this 
non monotonic behaviour is that the LLL regime of C2 
is reached at smaller values of Tlll with respect to the 
case of the vortex core (see Fig. . The results for the 
elastic shear modulus are relevant for the interpretation 
of the experiments of Ref. [3(| on the Tkachenko oscilla- 
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tions in harmonically trapped condensates. Indeed, one 
can average the predicted bulk values for Ci over the TF 
profile of the trapped condensate, using a local density 
approximation (see for example Ref. |29j). This proce- 
dure allows to estimate the Tkachenko frequencies also 
in the intermediate regime between the TF and the LLL 
limit. In practice, for the highest rotation rate realized 
in the experiments, Q/uj± = 0.99, where the gas is in a 
2D regime and co± is the radial trapping frequency, we 
find that C2 is lowered by only about 30% with respect 
to the TF value (n)/8. This means that the experimental 
data are still far from the LLL limit of the Tkachenko fre- 
quency, where the shear modulus would depend linearly 
on g 1 1 11 ] - Furthermore, the predicted reduction of C2 is 
not sufficient to explain the sizable discrepancy between 
the experimental measurements and the TF prediction 
|29| . Possible explanations of the remaining discrepancies 
are: (i) inadequacy of the local density approximation in 
the calculation of the Tkachenko frequencies, (ii) effects 
associated with the strong observed damping, (iii) de- 
viations from linearity in the experimental excitation of 
vortex modes, and (iv) occurrence of anharmonic effects 
in the trapping potential. 

In conclusion, we have investigated the bulk proper- 
ties of a rotating condensate in the presence of exact 
compensation between the centrifugal potential and the 
harmonic trapping. To this purpose, the mean field GP 
equation has been solved by introducing a gauge trans- 
formation which explicitly exploits the periodic structure 
of the density distribution and of the velocity field in the 
rotating frame. We have systematically investigated the 
role of interactions exploring the transition between the 
Thomas- Fermi and the LLL regime. We have found that, 
while the size of the vortex core exhibits a smooth mono- 
tonic behaviour as a function of the interaction parame- 
ter, the elastic shear modulus is characterized by an in- 
termediate maximum. The comparison with experiments 
reveals a good agreement as concerns the size of the vor- 
tex core and the static properties of the lattice, whereas 
we proved that the remaining discrepancies regarding the 
Tkachenko oscillations measured at the highest angular 
velocities cannot be attributed to LLL effects alone. 
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J. Dalibard, and N.H. Lindner for fruitful discussions. 
We would like to thank E. Cornell and V. Schweikhard for 
valuable comments and sharing their data. C.T. would 
like to warmly thank M. Paolini for introducing him to 
the finite element method. 
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